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ABSTRACT 

We present analytic and numerical models of the 'cluster wind' resulting from the 
multiple interactions of the winds ejected by the stars of a dense cluster of massive 
stars. Wc consider the case in which the distribution of stars (i. e., the number of stars 
per unit volume) within the cluster is spherically symmetric, has a power-law radial 
dependence, and drops discontinuously to zero at the outer radius of the cluster. We 
carry out comparisons between an analytic model (in which the stars are considered 
in terms of a spatially continuous injection of mass and energy) and 3D gasdynamic 
simulations (in which we include 100 stars with identical winds, located in 3D space 
by statistically sampling the stellar distribution function). From the analytic model, 
we find that for stellar distributions with steep enough radial dependencies the cluster 
wind flow develops a very high central density and a non-zero central velocity, and 
for steeper dependencies it becomes fully supersonic throughout the volume of the 
cluster (these properties are partially reproduced by the 3D numerical simulations). 
Therefore, the wind solutions obtained for stratified clusters can differ dramatically 
from the case of a homogeneous stellar distribution (which produces a cluster wind 
with zero central velocity, and a fully subsonic flow within the cluster radius) . Finally, 
from our numerical simulations we compute predictions of X-ray emission maps and 
luminosities, which can be directly compared with observations of cluster wind flows. 
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1 INTRODUCTION 

Super star clusters (SSCs) are dense clusters of young mas- 
sive stars, first identified in NGC 1705 by Melnick et al. 
(1985) and in NGC 1569 by Arp & Sandage (1985). Re- 
cently, they have been observed in a wide range of star- 
forming galaxies, such as merging systems (NGC 4038/4039, 
Whitmore & Schweizer 1995), dwarf galaxies (Henize 2-10, 
Johnson et al. 2000), classical starbursts (M82, Gallagher 
& Smith 1999, Melo et al. 2005), as well as in our galaxy 
(Arches Cluster : Nagata et al. 1995; Cotera et al. 1996; 
Serabyn, Shape & Figer 1998) amongst many other systems 
(for a review see Whitmore 2000). 

These star clusters can contain hundreds or thousands 
of very young, energetic stars, and have stellar densities far 
greater than those seen in normal OB associations. The ages 
of most of these star clusters are around 1- 10 Myr, their 
radii typically in the range of --^ 1-10 pc, and their total 
cluster masses in the 10^ - 10^ M© range (Melo et al. 2005 
reported a mean mass per star cluster, of ~ 2 x 10^ Mq, 
for M82). The central stellar densities of SSCs reach up to 
~ 10^ Mq pc~^. However, we can find SSCs with older ages 
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and/or larger masses (Walcher et al. 2006 reports a cluster 
with 6 X 10^ Mq). 

Canto et al. (2000, hereafter Paper I) explored an ana- 
lytic model using the physical properties of the Arches Clus- 
ter. This cluster has ~ 100 massive stars within a ~ 0.2 - 
0.3 pc radius. Such SSC must have strong multiple stellar 
wind interactions, resulting in the formation of stellar wind 
shocks. In Paper 1, the cluster wind was modeled as a mass- 
loading process (Hartquist et al. 1986; Dyson 1992; Lizano 
et al. 1996), resulting in essentially the same model as the 
one presented by Chevalier & Clegg (1985, who studied the 
flow resulting from a continuous temporal and spatial dis- 
tribution of supernovae explosions). 

Both in the stationary solution for spherically symmet- 
ric winds (Chevalier & Clegg 1985, Canto et al. 2000) and in 
the numerical calculations of Raga et al. (2001) the stellar 
distribution (within the cluster) was assumed to be homoge- 
neous. Also, these models are adiabatic (or, more precisely, 
non-radiative) solutions, which are appropiate for SSCs with 
low to intermediate mass and/or terminal velocity of SSCs 
(masses around 10'* - 10^ Mq and terminal velocities of 
~ 1000 km/s). For more massive stellar clusters, lower stel- 
lar wind terminal velocities or higher metallicities, radiative 
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losses within the cluster wind may become important (see 
Silich et al. 2004). 

In the present paper, we explore a non-radiative 
analytic model for cluster winds, considering a non- 
homogeneous stellar distribution (within the cluster). In 
particular, we study the case in which the stars have a 
n{R) oc R" power-law distribution (where n is the number 
of stars per unit volume as a function of the spherical radius 
R), with —3 < a ^ 0. We also compute 3D gasdynamical 
simulations for stellar distribution functions with different 
values of a, and compare the properties of the computed 
flows with the analytic model. This work is a natural exten- 
sion of the model presented in Paper I. 

We note that Matvienko & Shchekinov (2005) have pre- 
sented a study of mass loaded flows with sources and sinks 
with power law spatial distributions. The analytic model 
that we discuss in the present paper is based on a similar 
set of equations, but describes a cluster wind, rather than 
a mass loaded stellar wind (as was studied by Matvienko & 
Shchekinov 2005) . Also relevant in the context of the present 
work are the papers of Raga et al. (2001) and Rockefeller 
et al. (2005), who carried out 3D numerical simulations of 
winds from stratified clusters. 

The paper is organized as follows. In section 2, we 
present the analytical solution. In section 3, we describe 
the numerical simulations and compare the results with 
the analytic model. Predictions of the X-ray emission from 
the simulated flows are described in section 4. Finally, we 
summarize our results in section 5. 



2 THE ANALYTIC MODEL 

We consider A'' identical stars in a spherical cluster with 
an outer radius Re- The stars have a spatial distribution 
(number of stars per unit volume) of the form, 

(3 + Q)iV 



n{R) = k^R° 
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(1) 



where R is the spherical radius, a and fee are constants. 
In the second equality, the constant has been computed 
using the normalization condition 



A-KR^n{R) dR = N . 



(2) 



Every star has an identical wind with mass and energy 
depositon rates and Ew, respectively, and terminal 



velocity Vw For highly supersonic wind, Eu 



The stellar winds are thermalized at shocks (produced by 
interactions between the multiple winds), resulting in the 
production of a hot intercluster gas. This gas has a large 
central overpressure, which eventually produce a stationary 
'cluster wind' flow. 

For an adiabatic, spherically symmetric cluster wind 
(and negleting the gravity due to the stellar distribution), 
the mass, momentum and energy equations are : 
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p and V are the mass density and velocity of the wind re- 
spectively, R is the radial coordinate. 
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and his the speciflc enthalpy, P is the mean gas pressure 
and 7 is the specific heat ratio. 



Equation (|3]) implies that 
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where Mc = NMw is the total mass loss rate from the 
cluster in steady state. This equation shows that the mass 
flux profile strongly depends on the a exponent. 
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The adiabatic sound speed Cs is given by, 
fcT 
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where T is the gas temperature, p is the mean mass per 
particle and k is the Boltzmann constant. Integrating equa- 
tion ((5]) and using equations (|6} and ([8]) we find 



c2 = l- 
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Now, combining equations Q, ((Tjl and (|9]) we obtain 



(10) 



1 — cv dv dr 
a + bv^ ' 

where, a — 1 + a, b — (l + a+[5 + ia]7)/(7 — 1), and 
c = (7 + 1) /(7 — 1). In terms of the dimensionless variables, 
V — V/Vw r = R/Rc , this equation admits solutions, 

(^a + bv^^'' = Ar v = {—a/b)^^"^ = const. (11) 

In equation (|lip . p = 1/a, q = ~{b + a c) / (2a b) and A is an 
integration constant. 

Outside the cluster (i. e., for R > Rc) n{R) — 0, and 
the mass and momentum conservation equations have the 
form, 

V— - - — 

^ dR^~ dR' 



Combining equations (|12|l and (|13|l . we obtain 
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where, B is a constant. 

The velocity of the fiow, (wi) at the outer boundary of the 
cluster (r = R/Rc — 1) follows from equation (|14|l . 
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The value of the integration constants A and B are de- 
termined by the boundary conditions at r = 1 (7? = Rc) 
and r — > 00. They, of course, depend also on the value of a. 
There are three different regimes : 
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a. a > — 1 : 

In this regime, p > and g < (also both a and b >0). 
From equation we can see that very close to the centre 
of the cluster, 



(16) 



The wind velocity inside the cluster is subsonic, the veloc- 
ity at the centre of the cluster is zero and the density is 
finite. The gas velocity increases towards the outer cluster 
boundary, and the flow at this boundary {r — 1) follows 
from equation ( llip . 

A^v'-la + bviy. (17) 
The left-hand side of equation (|15p has a maximum value of 

1/2 / „ \ 1/(7-1) 

(18) 
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for vi = [(7— l)/(7+ l)]^''^- The requirement that the pres- 
sure at infinity must go to zero implies that the flow has to 
adopt the critical solution, for which 
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We can then use equations (O, ((8j, (O and (fT6|l to 
obtain the properties of the cluster wind at the centre of the 
cluster (see Paper I): 
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< a < -1 



In this interval, p<0, q>0, a<0 and b > 0. Thus 
the velocity inside the cluster is subsonic and has a non-zero 
value at the cluster centre (r = 0). From equation pip 
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On the other hand, equation ((71 indicates that the mass 
flux of the wind is inversely proportional to a positive power 
of r. This implies that the central gas density tend to infinite. 
From equations (|8]), (|9]) and (|24p the central temperature is, 



To = 



7 - 1 M 
27 k 



(-1) 



V 



(25) 



C. Oicr — 3 ^ ^ OLmin • 

In this regime, p>0, g<0, a<0 and 6 > 0, the veloc- 
ity of the wind remains constant (with a supersonic value) 
within the cluster. Its magnitud is given by equation (|24p . 
The temperature is also uniform inside the cluster and is 
given by equation (|25|l. The density also goes to infinity at 
the centre of the cluster and decreases outwards. But the 



wind escapes from the cluster surface supersonically, and 
accelerates outwards, until it reaches the terminal velocity. 

We do not consider power-law stellar distributions with 
a < Qci- {= —3) because they have an infinite number of 
stars (resulting from the strong divergence of the distribu- 
tion function in the cluster centre). 



3 THE NUMERICAL SIMULATIONS 
3.1 Numerical setup 

In order to illustrate the analytic star cluster wind solutions 
we have computed 3D numerical simulations with the 
full, non-radiative gasdynamic equations. The simulations 
solve a multiple stellar wind interaction problem with the 
3D, adaptive grid "yguazii-a" code, which is described 
in detail by Raga et al. (2000, 2002). The simulations 
were computed on a five-level binary adaptive grid with a 
maximum resolution of 0.1172 pc (corresponding to 256'^ 
grid points at the maximum grid resolution) in a computa- 
tional domain of 30 pc (along each of the 3 coordinate axes). 

In all runs, we assumed that the computational do- 
main was initially filled by a homogeneous, stationary am- 
bient medium with temperature Te„„ = 500 K and den- 



sity He 



= 0.1 cm" 



The stellar winds are imposed in 



spheres (centred at the stellar positions, see below) of radius 
Rw = 2.2 X 10^* cm, corresponding to 6 pixels at the maxi- 
mum resolution of the adaptive grid. Within these spheres, 
we impose (at all times) a Tu, = 15000 K temperature, and 
an outwardly directed Vw = 1000 km s~^ velocity. The den- 
sity within the spheres has an law (where r is the radial 
coordinate measured outwards from the stellar position), 
scaled so that the mass loss rate is Mn> = 10~^ Mq yr~^ 
for each star. We then place 100 such stellar wind sources 
within a spherical cluster of outer radius Tic = 10 pc, centred 
in the computational domain. 

We have computed four simulations of clusters with 
power-law stellar distributions (see equation [TJ with a — 
0, —0.5, —2.0 and —2.5. The stellar positions are chosen by 
statistically sampling the distribution functions in the way 
described in section [331 



3.2 Sampling the stellar position distribution 
function 

In order to produce models that can be compared with the 
analytic solutions, we place stars within the cluster with po- 
sitions obtained by randomly sampling the power-law dis- 
tribution function given by equation ((TJ. This is done as 
follows. 

We first note that f{R)dR = C A-KR^n{R)dR (with 
n{R) given by equation [T} is the fraction of the stars which 
have radial positions between R and R -f dR. Using the nor- 
malization condition " f{R)dR — 1, we obtain 

m) = (j^) R'^" ■ (26) 



This normalization is of course valid only for a > — 3 since 
otherwise a divergence occurs at the lower limit of the nor- 
malization integral. 
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With a random number generator we then choose a 
number r) which is uniformly distributed in the interval [0, 1] . 
This variable is statistically related to the radius R through 
the relation : 

R p 

f{R)dR= / d-q. (27) 
Jo 

From equations (|26I27|) we obtain 

i?cr?'/^'+"^ (28) 



R 



from which we can sample J? as a function of the random 
number rj. 

Once we have chosen the radial coordinate for each of 
the A'^ stars (used in the simulation) by sampling the radial 
distribution function (as described above), we asign random 
directions to the vector position of each of the stars. In this 
way, we obtain statistical samplings of the stellar position 
functions. 

In practice, we have to modify the obtained stellar dis- 
tributions because we impose the stellar wind conditions in 
spheres of a finite radius (see section 2.1). Whenever 
we obtain pairs of stellar positions resulting in overlapping 
"stellar wind spheres", we elliminate one of the two stars. 
This leads to an undersampling of the desired distribution 
function in the central, high stellar density regions of the 
generated clusters. Because of this, the comparison between 
the numerical simulations and the analytic cluster wind so- 
lutions is only meaningful away from the central region of 
the cluster. 

We should note that the 100 stars of each model were 
chosen with the same set of random numbers, so that the 
directions from the cluster centre to each of the stars are the 
same in all models. However, the physical radii correspond- 
ing to the "radial random numbers" differ for each model, 
as the "conversion" from random number {rj) to physical 
radius R (see equation [28} depends on the value of a. 

3.3 Model results 

We computed four numerical simulations of clusters with 
stellar distributions with different a values. All models are 
computed for clusters with 100 stars and a 10 pc outer ra- 
dius. The stars are identical, with a mass deposition rate of 
10~^ Mq yr~^ and a wind velocity of 1000 km/s. 

The only difference between the models are the stellar 
density distributions. We have chosen a values (see equa- 
tion [l]) covering the different regimes of the cluster wind 
flow (see section [2]). 

We have integrated forward in time the four models un- 
til stationary flows are obtained. From the stationary flow 
configurations, we then compute (through appropriate inter- 
polations in the cartesian adaptive grid) the radially depen- 
dent flow density, velocity and temperature averaged over 
spherical concentric surfaces Sr — A-kR? : 



P'^i^) = ^ / P ^^^^ 



Va{R) = 
Ta{R) 



4npa{R)jg 

1 

ATipa{R)jg 



pVR sin 9 d9d(j> , 
pT sin6l( 



(29) 
(30) 
(31) 




Figure 1. The stellar distribution of xy-plane and xz-plane for 
the a = model. The solid line represents the outer boundary of 
the cluster. 



where 9 and (j> a-r^ the polar and azimuthal angles, respec- 
tively, and p is the flow density, T the temperature and vr 
the radial velocity (obtained by projecting the three carte- 
sian velocity components resulting from the numerical inte- 
gration onto the direction normal to the spherical surface). 
That is Vr — {xvx + yvy -\- zVz)/R- We then compare this 
spherically averaged flow with the results from our analytic 
cluster wind model. 



3.3.1 The -1 case 

For this a range we ran two models: a cluster with a 
homogeneous (i.e., a = 0) stellar density distribution, and 
a cluster with an a = —0.5 distribution. 

The stellar distribution used for the a = Q model (ob- 
tained by sampling the a = distribution function, see equa- 
tion [21] and section I3.2|l is shown in Figure [T] The average 
flow variables as a function of spherical radius (obtained 
by averaging over concentric spherical shells, see above) 
are compared with the solution obtained from the analytic 
model in Figure [2] From this Figure, it is clear that the av- 
erage density, velocity and temperature (obtained from the 
numerical simulation) agree very well with the analytical 
model, except close to the centre of the cluster. The devi- 
ations in the cluster centre are a direct result of the small 
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^ variables computed from the flow that results from the nu- 

'° '5 merical simulation, except for the inner region of the cluster. 

Radius [pc] 



Figure 2. The spherically averaged flow obtained from the a = 
model. The density (top), radial velocity (centre) and temper- 
autre (bottom) obtained from the numerical simulations (dashed 
lines) and the analytic model (solid lines) are shown as a function 
of spherical radius R. 

number of stars that is present within the inner ~ 3 pc of 
the cluster (see Figure [TJ. 

Figure [3] illustrates the real complexity of the flow that 
is obtained from the numerical simulations. This Figure 
shows a 3D rendition of the 10 pc radius star cluster, with 
100 stars and a = 0. The flow has a sponge-like morphol- 
ogy, with low density stellar wind cavities immersed in the 
denser 'cluster wind' flow (composed of shocked stellar wind 
material). The relatively monotonic average flow (shown in 
Figure [2} is obtained by averaging over this complex flow 
structure. 

The stellar distribution used for the a = —0.5 model 
(obtained by sampling the a = —0.5 distribution func- 
tion, see equation 1261 and section is shown in Figure 
The average flow variables as a function of spherical radius 
(obtained by averaging over spherical surfaces, see above) 
are compared with the solution obtained from the analytic 
model in Figure [5] Again, a reasonably good agreement is 
obtained between the analytic model and the average flow 



3.3.2 The Umin < a < ac case 

Figure |6] shows the stellar positions that results from sam- 
pling the Q = — 2 distribution function. The stars are highly 
concentrated around the centre of the cluster. 

Figure [7] shows the angularly averaged velocity, density 
and temperature as a function of the radial coordinate. The 
analytical solution has a non-zero, subsonic velocity and an 
infinite density at the central position of the cluster. 

Substantial differences between the numerical and 
analytical solutions are found for radii smaller than ~ 4 pc. 
This is a direct result of the undersampling of the stellar 
distribution function which occurs as a result of the 
"proximity criterion" (described in 13. 2|) applied for placing 
the stars in the computational grid. For larger radii, a 
reasonable agreement between the analytical and numerical 
results is obtained. 



3.3.3 The Ucr < a < ctmin case 

Figure |8] shows the stellar positions that results from sam- 
pling the a = —2.5 distribution function. The analytical 
solution and the angularly average flow variables (obtained 
from the numerical simulation) are shown in Figure |9] For 
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Figure 3. 3D rendition of tlie model with 100 stars distributed homogeneously (a = 0) inside of a sphere of 10 pc in radius. In logarithmic 
grey-scale (colour-scale in the online version) we present 5 isosurfaces of density. Depicted by arrows we overlaid the velocity field, the 
largest arrow corresponds to a magnitude of 10^ km . 



this model, the analytic solution has a supersonic, outwards 
velocity in the inner region of the cluster. Again we ob- 
tain substantial deviations between the analytic and numer- 
ical solutions in the central region of the cluster, and better 
agreement for larger radii. 



4 PREDICTED X-RAY EMISSION 

We have taken the density and temperature stratifications 
for the steady cluster wind flow configurations, and used 
them to compute the X-ray emission. We have done this 
by computing the emission coefficient in the 0.3 2 keV 
photon energy range using the CHIANtJ] atomic data base 
and software (see Dere et al. 2001 and references therein). 

^ The CHIANTI database and associated IDL pro- 
cedures, now distributed as version 3.0, are freely 
available at the following addresses on the World 
Wide Web: |http:/ /wwwso lar.nrl. na vy.mil/chianti. html[ 



For this calculation, it is assumed that the ionisation state of 
the gas corresponds to coronal ionisation equilibrium, and 
that the emission is in the low density regime (i. e., that 
the emission coefficient is proportional to the square of the 
density) . 

We then integrate this emission coefficient along lines of 
sight, which are assumed to be parallel to the j/-axis of the 
computational grid. The X-ray maps computed in this way 
for the models with a = 0, —0.5, —2 and —2.5 are shown in 
Fig. 10. 

From these maps, it is clear that the X-ray emission is 
highly structured in all models. Actually, some of the struc- 
tures are seen in all four models, as the stellar position have 
been chosen with the same set of random numbers (but not 
the same radial positions, see the discussion at the end of 
§3.2). 

In Figure 11, we show the emission along a line of sight 



http: //www. arcetri. astro. it/science/chianti/chianti.htmT| an d 
http : / / www . damtp . cam .ac.uk/user /astro / chianti / chianti .html] 
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Figure 10. X-ray images obtainded by integrating the emission coefficient along lines of sight (parallel to z-axis) and over energies 
ranging from 0.3 to 2 keV, for the four models described in the text. The images arc depicted witii gray scale (with the fluxes normahzcd 
to the maximum value) given by the bar to the top of the right-top plot). The scales along the x- and y-axes are given in units of parsecs, 
with the origin coinciding with the barycenter of the stellar distribution of the cluster. 



placed at x, j/=0 four X-ray maps of Figure 10. These cuts 
show that the a = and a — —0.5 models have a structured 
but more or less constant emission within a radius of ^ 6 pc 
(i. e., the inner 60 % of the cluster radius), with wings of 
lower emission extending to beyond the outer cluster bound- 
ary. The models with steeper stellar distribution functions 
(i. e., with more negative values of a) show more compact 
emission structures. 

Finally, by integrating the X-ray emission over the 
whole emitting volume, we compute luminosities of 0.32, 
0.32, 0.26 and 0.20 L© for the a = 0, -0.5, -2 and -2.5, 
respectively. Therefore, we find that the while the predicted 
maps differ quite appreciably between the different models, 
the X-ray luminosity does not depend strongly on the value 
of a. 

We should note that because the emission coefficient 
is proportional to the square of the density, one can use 
equation (21) to deduce the scaling of the X-ray luminosity 



to other model parameters. For example, a supcrclustor with 
~ 20000 stars would have X-ray luminosities of ~ 8 — > 12 x 
10^ Lq (depending on the value of a, see above). 



5 CONCLUSIONS 

In this paper we have extended the analytic cluster wind 
model of Paper I to the case of a non-uniform stellar distri- 
bution. In particular, we have studied the case of a radially 
dependent, n{R) — kcR" power-law distribution. 

Of particular interest is the a = —2 distribution, which 
corresponds to the stratification of a singular, isothermal, 
self-gravitating sphere. Such a structure is of interest for 
modelling the wind from a gravitationally bound stellar clus- 
ter. Power-law stellar distributions with other values of a do 
not have a clear physical justification, but can be considered 
as a parametrization of stellar distributions with different 
degrees of central condensation. 
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Figure 5. Same as Figure[2l but for the a = —0.5 model. 



We find that for shallow distributions, with — 1 < a < 
0, the cluster wind has zero velocity in the cluster centre. 
However, for more negative a values the cluster wind has a 
non-zero, outwards directed velocity (subsonic for amin < 
a < —1 and supersonic for — 3 < a < amin) in the centre of 
the cluster. The solutions with a non-zero central velocity 
have an infinite central density for the cluster wind. 

We have then compared the analytic cluster wind solu- 
tions with 3D numerical simulations. For carrying out the 
simulations, we consider the winds from 100 stars, with a 
spatial distribution obtained by statistically sampling the 
appropriate stellar distribution function. We then carry out 
angular averages of the computed flow variables, and com- 
pare the radial dependence of these averages with the pre- 
dictions obtained from the analytic model. 

For different values of a, we obtain a good agreement 
between the analytic and numerical predictions in the outer 
regions of the cluster. However, the analytic and numerical 
solutions have large differences in the central region of the 
cluster. These differences are a direct result of the fact that 
only a small number of stars are present in this spatially 
reduced region, and therefore the continuous mass and en- 



ergy source distribution assumed in the analytic model is 
inappropriate for describing the real cluster wind flow. 

Finally, we have obtained predictions of the X-ray emis- 
sion from our simulated cluster wind flows. Our models have 
quite low ~ 0.2 0.3 L© X-ray luminosities, as a result 
of the fact that only 100 stars were included (due to the 
constraints imposed by the numerical resolution of the sim- 
ulations). However, the scaling laws of the analytic model 
imply that models with 20000 stars (i. e., with the number 
of stars of a large supercluster) will have X-ray luminosities 
of ~ lO** Lq (see section 4). At distances of ~ 4 Mpc (i. e., 
the distance to M 82), this would produce an X-ray flux of 
2 X 10"" erg s"^ cm"^. Such a flux is only 2 % of the X- 
ray flux observed for the Arches cluster close to the galactic 
centre (see Yusef-Zadeh et al. 2002), but might be within 
range for possible future observations. 

An important remaining problem is that Ha emission 
has been observed from many superclusters. Our models 
have gas temperatures in excess of lO'^ K, and therefore do 
not produce such an emission. Therefore, the observed Ha 
emission must be coming from another component, which 
could be the remnants of dense interstellar clouds which were 
present in the region when the cluster was formed. The pres- 
ence of such dense, lower temperature structures within the 
cluster wind flow could be explored in the future with nu- 
merical simulations similar to the ones which we have shown 
in the present paper, but that need to include the radiative 
cooling. 
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Figure 7. Same as Figure |2] but for the a = —2 model. 
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